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Abstract— Considering the characteristics of the x-ray mi- 
croscope system being developed at SUNY at Buffalo and the 
limitations of available cone-beam reconstruction algorithms, 
a general cone-beam reconstruction algorithm and its several 
special versions are proposed and validated by simulation. Our 
cone-beam algorithm allows various gowning lod, handles recon- 
struction of rod-shaped specimens which are common in practice, 
and facilitates near real-time reconstruction by providing the 
same computational efficiency and parallelism as Fddkamp's 
algorithm [1L 

Although our cone-beam algorithm is not exact, It consistently 
gives satisfactory reconstructed images. Furthermore, it has sev- 
eral nice properties if the scanning locos meets some conditions, 
First, reconstruction within a mldplane is exact using a planar 
scanning locus. Second, the vertical integral of a reconstructed 
image is equal to that of the actual image. Third, reconstruction is 
exact if an actual image is independent of rotation axis coordinate 
z. Also, the general algorithm can uniformize and reduce z-axis 
artifacts, if a helix-like scanning locos is used. 



I. INTRODUCTION 

AN X-RAY shadow projection microscope system using 
a scannable x-ray point source is under development at 
AMIL-ARTS (SUNY at Buffalo) for acquiring three-dimen- 
sional micro-images, which are very useful in the study of 
biological and synthetic materials [2H9] (Fig. 1). The x-ray 
microscope system uses an electron optical column to produce 
a small, intense, and steerable electron spot, which in turn 
creates an x-ray spot of sufficient brightness to illuminate a 
specimen to resolve the structures of interest. The resolution 
of this x-ray microscope system is limited by the size of the x- 
ray spot due to penumbra! blurring, which can be minimized 
by reducing the lateral spread of the electron beam in the 
target material [3]. The specimen is on the side of the target 
away from the electron beam. This permits the electron beam 
to be in vacuum while the specimen can sit in air or in an 
environmental chamber. To obtain the needed magnification, 
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Fig. 1. The schematic diagram showing the x-ray shadow projection 
microscope being developed at SUNY at Buffalo. 

the specimen must be placed close to the source, in this case, in 
the range of 0.1-1 mm. An x-ray projection image is obtained 
by placing a phosphor screen as the x-ray scintillator at a 
large distance from the specimen. The larger the distance, the 
higher the magnification, but the lower the image intensity. 
A two-dimensional slow-scanned CCD array (Kodak 1400) is 
used to grab the projection image formed on the screen. For 
microtomography, a series of such images is acquired at dif- 
ferent orientations, and then processed for reconstruction. This 
approach is commonly called cone-beam reconstruction, as x- 
rays irradiating to a sphere-like specimen form a cone (Fig. 2). 

Because none of the currently available cone-beam re- 
construction algorithms are fully satisfactory for our x-ray 
microscope system, we will propose in this paper a general 
cone-beam algorithm which can be considered as a gener- 
alization of Feldkamp's algorithm [1]. Currently available 
cone-beam reconstruction algorithms are briefly reviewed in 
Section H. A heuristic mathematica] derivation of the proposed 
general cone-beam algorithm is given in Section m. In Section 
IV, several interesting special versions of our cone-beam 
algorithm are studied. Analytic properties of our cone-beam 
algorithm are discussed in Section V. 

II. Review of Cone-Beam Algorithms 

The cone-beam reconstruction of x-ray images has been 
studied for more than thirty years. The algorithms can be eval- 
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Fig. 2. Cone-beam reconstruction arrangement 

uated from various viewpoints, such as theoretical precision, 
computational efficiency, scanning geometry, and reconstruc- 
tion region. A recent review was written by Smith [10]. 

Existing cone-beam algorithms can be categorized into two 
broad classes: analytic and algebraic [11], [12]. Algebraic 
algorithms demand more computing resources than analytic 
ones. Because a high reconstruction speed is very important 
in practice, analytic algorithms are always preferred. 

Kirillov [13] developed an analytic formula for recon- 
structing a complex valued function in an n-dimensional 
complex space using complex valued cone-beam projection 
data. Kirillov assumed that the vertex describes an unbounded 
curve and showed that a sufficient condition for reconstruction 
in the Schwarts space is that the curve intersects almost every 
hyperplane. Tuy [14] modified Kirillov' s approach and gave a 
reconstruction formula for functions of compact support when 
the curve is bounded, provided that almost every hyperplane 
meeting the support of the unknown function meets the 
curve at some point Tuy derived a cone-beam reconstruction 
formula for two intersecting source circles, which requires a 
gradient to be computed at each vertex point. Hamaker et 
aL [15] discussed an approximate approach for cone-beam 
reconstruction when there are a finite number of source points 
available. These results are theoretically oriented and not 
directly applicable in practice. 

Nalcioglu and Cho [16] and Denton et al [17] presented 
convolution-backprojection methods requiring that source po- 
sitions encompass a sphere. The analytic results obtained by K. 
T. Smith [18], [19], Imiya [20], and Peyrin [21] are also based 
on the assumption of a sphere of 'Vertex points.** However, it 
is very difficult to implement such a geometry, particularly in 
our x-ray microscope system. 

Herman [22], and Lewitt and McKay [23], studied the 
unmodified fan-beam algorithm for cone-beam reconstruction. 



When the cone angle is small as in the scanner developed 
at the Mayo Clinic [24], the convolution-backprojection pro- 
cedure works quite well. Grangeat's algorithm [25], [26] is 
also in the form of convolutiori-backprojection, where the 
planar integral values are estimated for three-dimensional 
Randon inversion. Then, Sire, Grangeat, and Rizo et al [27] 
developed a software package named RADON. Feldkamp 
et al [1] proposed a practical cone-beam algorithm, which 
is an extension of the conventional equispatial fan-beam 
reconstruction algorithm, where two-dimensional projection 
data from different angles are filtered and backprojected to 
voxels along projection beams. The final value of each voxel 
is the sum of contributions from all tilted fan-beams passing 
through the voxel. Jaszczak's group and Gullberg's group also 
studied approximate filtering and backprojection cone-beam 
reconstruction [28H31] 

Guliberg etal [28], [29] and Wang etal, [32] independently 
studied Feldkamp's algorithm for planar noncircular scanning 
loci in different ways. Wang's hueristically derived formula 
is the same as Feldkamp's, except that the source-to-origin 
distance depends on the rotation angle, which is a special case 
of our cone-beam algorithm presented in this paper. Gullberg's 
formula contains a derivative of scanning locus resulted from 
the Jacobian transform, while Wang's formula does not contain 
such a derivative. 

A thorough theoretical analysis on cone-beam reconstruc- 
tion was done by Smith [10], [33]-[36]. He contributed 
the necessary and sufficient conditions for exact cone-beam 
reconstruction and the corresponding reconstruction methods. 
His algorithms are generally used in "complete** geometries 
and are less efficient than Feldkamp's algorithm, which is 
applied to data collected from a circle, a typical case of incom- 
plete scanning geometries. Although his most computationally 
efficient algorithm is as efficient as Feldkamp's algorithm, it 
is based on the following assumption of imaging geometry: 
Considering any line that contains a vertex point (a position of 
the x-ray source) and a reconstruction point, there is an integer 
M (which remains constant for this line) such that almost every' 
plane that contains this line intersects the geometry exactly 
M times. Known practical loci cannot meet this requirement. 
Actually, the only known curve mat satisfies this condition is 
an infinitely long straight line extended in both directions. 

Among various approaches, Feldkamp's heuristically de- 
rived algorithm is the most efficient and most suitable to 
parallel processing. For cone-beam x-ray microtomography, 
the main drawbacks of Feldkamp's algorithm, shortcomings 
of most of other algorithms as well, are that a specimen must 
be contained in a sphere-like reconstruction region and that 
the x-ray source must be moved along a circle in the spec- 
imen coordinate system. In addition, Feldkamp's algorithm 
suffers from a larger inaccuracy in reconstructing off-midplane 
structures. 



III. Heuristic Derivation of 
General Cone-Beam Algorithm 

Our cone-beam algorithm is derived by modifying the 
convolution and backprojection algorithm for equispatial fan- 
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beam data. The locus of the x-ray source is described in 
a cylindrical coordinate system (p{0) 1 h{fl) 1 0) f where 0 € 
[0,B] is the rotation angle, p(0) describes the horizontal 
distance between the source and the origin of the coordinate 
system, h{0) denotes the vertical distance from the source 
to the reference polar coordinate plane. The cylindrical co- 
ordinate system can be converted to the Cartesian coordinate 
system using 



X = (>(0)CO80 

* = h(0). 



0) 



It is assumed that the rotation axis is the z axis and the interval 
length of p is not less than 2tt, that is, B > 2*, which means 
that a specimen is illuminated from all horizontal directions. 

In our cone-beam reconstruction, there are mainly two 
scawiing modes, planar and helix-like. A typical case of planar 
scanning loci is a circle, which is used in the conventional 
cone-beam reconstruction. A helix-like scanning locus should 
be used to deal with rod-like specimens. Without lose of 
generality, a locus turn of the x-ray source can be described 
by the following parametric equations: 

[h = h(0) < 2 > 

where 0 e [0,2irJ. In the first scanning mode, k(0) = 0,0 e 
[0,2*-). In the second mode, h{0) € [z - + 4*],/? e 
[0,27r], where z is the longitudinal axis coordinate of a 
specimen transverse slice, and Ap is the pitch of a helix-like 
scanning locus. We define the midplane as the reference polar 
plane in the first scanning mode and the concerned horizontal 
plane with a longitudinal axis coordinate z in the second mode. 

The fan-beam reconstruction algorithm [12] can be written 
in the coordinate system (ar, y) as 

P) >J? 5 ° ,<W (3) 



where (*, s) denotes the rotated coordinate system described by 
* = xcos£ + y sin/3, s = -*sin/? + ycos/? 

R p(l>) represents the projection data, p describes a linear 
detector array intersecting the origin of the reconstruction 
coordinate system, D so denotes the distance between the 
x-ray source and the origin of the reconstruction coordinate 
system, and /(•) is the reconstruction filter defined as 



■c 



The specimen is assumed to be centered at the origin and 
along the z-axis (if it is rod-like). 

In the three-dimensional coordinate system (x,y,z) f the 
cone-beam reconstruction is performed by tilting the fan out 
of the midplane, thus changing the coordinate system (t, s) 
A new coordinate system is defined to represent the 
location of a reconstructed point with respect to the tilted fan. 




Fig. 3. Schematic diagram illustrating various coordinate systems and 
variables. 



Because of the geometry of the tilted fan, bom the source-to- 
origin distance and the angular differential are changed (Fig. 
3). the new source-to-origin distance p(0) can be determined 



P 2 (0) = P 2 (0) + < 2 



where C is the height of the fan above the center of the 
midplane as given below 



.- Mm 
tip)-*- 



(5) 



The local coordinate system associated with the 

angle 0 is used whose z coordinate is computed according to 

m = z-H0). 

Because the fan has been tilted, the equivalent angular incre- 
ment d0 for the tilted fan must be found such that 

d0p{0) = d0p(0). 

Then, the equivalent increment of angular rotation d0 is 

d0p{0) 



dp = 



(6) 



Substituting p{0) and d0 for D so and d0 respectively in 
(3) and denoting two^umensional projection data as R 0 (p } <), 



T 
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we obtain the incremental contribution dg(x,y>z) to the 
reconstructed voxel g(x, y , z) as 

or 



- S " v 'P 2 (/ J ) + P a + C 2 



Since 



*0 rtfl 



(9) 



by summing up all incremental contribution, we obtain the 
following formula for g(x,y,z): 



•/( 



(10) 



Equation (10) can be referred to as the generalized Feidkamp's 
cone-beam reconstruction formula [37], [38]. The voxels that 
are irradiated from all horizontal directions form the recon- 
struction region where reliable reconstruction can be achieved. 

The advantages of our generalized cone-beam algorithm are 
summarized as follows. First, this algorithm allows various 
scanning loci, giving much more freedom in mechanical and 
electrical implementation of an x-ray microscope system. 
For example, it allows a polygonal scanning locus, which 
eliminates the need of continuously rotating a specimen. 
Second, it applies to not only sphere-like specimen recon- 
struction but also rod-shaped specimen reconstruction, which 
is very common in practical applications. Third, it provides the 
same computational efficiency and parallelism as Feidkamp's 
algorithm. Finally, 2-axis artifacts can be uniformized and 
reduced by our generalized algorithm, if a helix-like scanning 
locus is used. These points will be illustrated in the following 
sections. 

IV. Numerical Simulation of 
General Cone-Beam Algorithm 

In this section, the digital implementation of our cone-beam 
reconstruction algorithm is explained. The simulation data 
and parameters are described. Then, several interesting cone- 
beam reconstruction algorithms are formulated and simulated 
as special cases of our cone-beam algorithm. 

A Digital Implementation of Our Cone-Beam Algorithm 

The key question in the digital implementation of our cone- 
beam algorithm is how to implement the reconstruction filter 
/(•) [12]. Let r denote the sampling interval of projection 



data. This means that in Fourier domain the energy outside 
the frequency interval (-W,W) is negligible where 



(ID 



Since the highest frequency of the projections is finite, we have 
R 0 (wX)\w\e j2irwt dw 

R 0 (wX)F(w)^ wt dw (12) 



where Rfl(wX) denotes the one-dimensional Fourier trans- 
form of Rp(p, 0 with respect to p. 



F(w) =\w\bw(w) 



and 



\w\<W 
otherwise. 



Therefore, f(t) can be determined as the inverse Fourier 
transform of F(w) as follows [12]: 

f 1 sin(7Tt/r) 1 , sin(7Tt./2T) x2 
,{t) ~ 2^ irt/r ~ ict/2r > ' CU) 

Since the projection data are measured with a sampling interval 
of t, f(t) needs only be known at discrete points with the same 
sampling interval. The samples, /(nr), of /(t) are given by 



/(nr) = < 0, n is even; 
l-^rfer, n is odd. 



(14) 



B. Simulation Data and Parameters 

A three-dimensional version of Shepp and Logan's head 
phantom {12] was used in the simulation. The head phantom 
parameters are listed in Table I. There are ten ellipsoids in 
the phantom. In Table I, xo>yoi26 specify the center of an 
ellipsoid, a y b,c are the x t y 3 z semi-axes respectively, a is 
the rotation angle of an ellipsoid (about z axis), and /i is the 
relative x-ray linear absorption coefficient. The effective x- 
ray absorption coefficient of a point is the sum of the relative 
parameters of ellipsoids containing that point Fig. 4 shows 
four slices of Shepp and Logan's three-dimensional head 
phantom. Two are horizontal slices with z = -0.25 (upper- 
left) and z = 0.625 (upper-right). Two are vertical slices with 
y = -0.105 Qower-left) and y = 0.1 (lower-right). Note 
that the real gray level was transformed according to a linear 
relationship for a better visualization. This transformation 
linearly maps the interval [0.95, 1.05] into 256 discrete gray 
levels. Truncation was done when necessary. 

In the numerical simulation, the detector plane is 2.2 by 
2.2 with 128 by 128 pixels and was so placed for the ease of 
computation that its center is at the origin of the reconstruction 
coordinate system and that its vertical axis is the z-axis of the 
coordinate system. One hundred projection images with a 3.6° 
angular interval were used for the reconstruction, except for the 
random cone-beam reconstruction where the angular position 
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TABLE I 

Parameters op Shepp and Logan's Threb-Dmensjonal 
Head Phantom Used en the Simulation 

No- ~xq to a K c » 



1 
2 
3 
4 
5 
6 
7 
8 
9 
10 



0.00 0.000 

0.00 0.000 

-0.22 0.000 

0.22 0.000 

0.00 0.350 

0.00 0.100 

-0.08 -0.650 

0.06 -0.650 

0.06 -0.105 

0.00 0.100 



0.000 
0.000 
-0.250 
-0.250 
-0.250 
-0.250 
-0.250 
-0.250 
0.625 
0.625 



0.6900 
0X624 
0.4100 
0.3100 
0.2100 
0.0460 
0.0460 
0.0460 



0.0560 



0.920 
0.874 
0.160 
0.110 
0.250 
0.046 
0.023 
0.023 
0.040 
0.056 



0.900 0 

0X80 0 

0.210 108 

0.220 72 

0.500 0 

0.046 0 

0.020 0 

0.020 90 

0.100 90 

0.100 0 



2.00 
-0.98 
-0.02 
-0.02 
0.02 
0.02 
0.01 
0.01 
0.02 
-0.02 



Fig. 5. Foot reconstructed slices obtained using Feldkamp'g cone-beam 
reconstruction algorithm when po = 3. corresponding to the four slices of 
the phantom shown previously. 






Rg. 4. Four slices of Shepp and Logan's three-dimensional head phantom 
= 7 0 - 25 ("PI** left). * = 0.625 (upper right), y = -0.105 (lower 
left), and y = 0.1 Cower right), respectively. 

is random. Feldkamp's cone-beam reconstruction algorithm is 
a special case of our algorithm with 



I P = Po 
I /i = 0. 



(15) 

Figure 5 shows the reconstructed images using Feldkamp's 
algorithm with Po = 3, corresponding to the four slices 
shown in Fig. 4. The averaged absolute errors between the 
original slices and the reconstructed ones were computed as 
3.5 (upper-left), 13.3 (upper-right), 13.2 (lower-left), and 13.3 
(lower-right) gray levels, respectively. 

C. Polygonal Cone-Beam Reconstruction Algorithm 

A polygonal cone-beam reconstruction algorithm is obtained 
[32] by 



/i = 0 



(16) 



where po is the minimum distance between the source and the 
origin, Ns is the number of sides of the scanning polygon, 
and [\J represents a floor function. The polygonal cone-beam 
reconstruction algorithm is used to reduce the number of 
recalibrations of the mechanical axis position. 



For tomographic imaging, a large number of two-dimens- 
ional projection images must be recorded at various ori- 
entations. The specimen therefore must be rotated, and a 
mechanical stage is installed for this purpose. Because the vari- 
ation of rotation axis position of the mechanical stage is one 
of the major error sources in microtomography, an elaborate 
feed-back system is required to minimize the error. However, 
experiments showed that the calibration of the rotation axis 
position is both time-consuming and complicated. On the 
other hand, the x-ray source of our microscope system can be 
precisely moved by steering the e-beam electromagnetically. 
With the polygonal algorithm, the number of the stage rotation 
can be greatly reduced. The straight line portions of the 
polygon are achieved by the x-ray source movement and the 
number of the polygon sides determines the number of the 
stage rotations. 

Figure 6 shows four reconstructed images obtained using 
the polygonal cone-beam reconstruction algorithm with po = 3 
and N s = 8. Hie averaged absolute errors between the original 
and reconstructed slices are almost the same as those obtained 
with Feldkamp's algorithm. We have proved that the polygonal 
algorithm will give exact reconstruction on the midplane, if the 
scanning polygon has an even number of sides. 

D. Multiple Cone-Beam Reconstruction Algorithm 

In this algorithm, Feldkamp's or the polygonal cone-beam 
algorithm is applied on different sets of cone-beam projection 
data captured with a set of planar scanning loci (Fig. 7) 
[39]. We only use the polygonal algorithm here. In general, 
if n x-ray sources are available, the interested reconstruction 
region can be horizontally divided equally into n segments 
and the central horizontal plane of each segment contains 
a polygonal scanning locus. Let us consider a twin cone- 
beam reconstruction algorithm. For a given slice with a 
vertical coordinate z, the closest planar scanning locus can 
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Fig. 7. The scanning locus in the multiple < 




Fig. 6. Four reconstructed slices obtained using the 
reconstruction algorithm with po — 3 tnd N$ = 8, 
slices of the phantom shown previously. 



corresponding to the four 



be described as 



P = 



(17) 



where po is a constant, JV 5 is the number of the scanning poly- 
gon sides, d is a vertical distance between the scanning planes, 
and sign(-) represents the sign function. This algorithm is used 
to improve the off-mid-plane reconstruction accuracy and to 
deal with rod-shaped specimens, which are commonly found 
in biological and material science. 

Two reasonable strategies for the estimation of voxel values 
are (1) to interpolate two voxel values obtained from the 
upper- and lower-cone-beam reconstructions; and (2) to select 
the voxel value obtained from the cone-beam reconstruction 
whose midplane is the nearest to the voxel. In the simulation, 
we used the second. Figure 8 shows the reconstructed images 
obtained using the twin polygonal cone-beam algorithm with 
po = 3, Ns = 8, and d — 1.25, respectively. Note that two 
midplanes used here are planes z = 0.625 and z ~ -0.625, 
while the midplane used for Feldkamp's algorithm and the 
polygonal algorithm is the plane z — 0. The lower two images 
show clearly horizontal dark stripes. The averaged absolute 
errors between the original and reconstructed slices are 5.3 
(upper-left), 2.4 (upper-right), 7.1 (lower-left) and 7.2 flower- 
right) gray levels respectively. This multiple algorithm reduces 
the reconstruction error of vertical slices, which is due to 
the fact that the reconstruction of off-midplane structures is 
improved. 

E. Helical Cone-Beam Reconstruction Algorithm 

In order to handle rod-shaped specimens, a helical cone- 
beam reconstruction algorithm can be designed. Contrast to 
Feldkamp's algorithm where the scanning locus of an x- 
ray source is a circle, the helical cone-beam reconstruction 
algorithm requires the translation of the x-ray source and the 



Fig. 8. 
cone-beam 
respectively, 
previously. 




slices obtained using 1 
algorithm with po = 3, N$ 
to the four slices of tJ 



i twin polygonal 
8. and d = 1.25. 
phantom shown 




Fig. 9. The scanning locus in the helical approach and a mid-plane. 

simultaneous rotation of the specimen, hence the motion locus 
of the x-ray source is a continuous helix in the specimen 
coordinate system. A helical turn (Fig. 9, right) is referred to as 
the portion of a helix generated by varying P continuously for 
2?r. For any horizontal slice of the specimen, it is assumed that 
such a helical turn can always be found that the turn is divided 
by the concerned slice into upper and lower halves. Projection 
images collected from this turn are used to reconstruct the 
specimen slice (Fig. 9.) A scanning turn in this scanning mode 
can be formulated as 



P = Po 



h 9 0 



(18) 



where po is a constant and h p is the pitch of the helix turn [40]. 

Compared with the multiple cone-beam algorithm, the heli- 
cal cone-beam algorithm results in a more uniform error dis- 
tribution along the z axis. Figure 10 shows four reconstructed 
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Kg. 11. The scanning locus in the broken-line 



approach and a mid-plane. 



Fig. 10. Four reconstructed slices obtained using the helical cone-beam 
n^nstroctton algorithm with p = 3 and h p I 1^5. con^o^to^e 
four sbces of the phantom shown previouslj. ""irsponoing to the 

images obtained using the helical cone-beam reconstniction 
algonto with po = 3 and hp = 1.25. The vertical distance 
K ~ l * 5 > wtween adjacent scanning planes in the multiple 
cone-beam algorithm is the same as the pitch of the helical 
scanning locus. It can be observed that the horizontal dark 
stripes shown in two lower images of Fig. 8 have been 
eliminated. The averaged absolute errors between the original 
and r^nstmcted slices are 4.3 (upper left), 3.4 (upper right), 
6^8 (lower left) and 6. 1 (lower right) gray levels, respectively 
These data show that the helical algorithm does uniformize 
and reduce the reconstruction errors between the original and 
reconstructed vertical slices. The error for the upper-right 
image (= 3.4) is higher than that (= 2.4) of the multiple cone- 
beam algorithm, because this slice is on one of the midplanes 
in the multiple cone-beam reconstruction algorithm. 

F. 3 D Broken line Cone-Beam Reconstruction Algorithm 

nie above helical cone-beam reconstruction utilizes a 
smooth and simultaneous rotation and translation of the 
specimen with respect to the incident x-ray cone. A 
three-dimensional broken-line cone-beam algorithm can be 
constructed from the combination of the polygonal cone-beam 
algorithm and the helical cone-beam algorithm. A turn of the 
scanning locus required by this hybrid algorithm is described 
by 



{ 



P = 



(19) 



coa(0-^ L £££ +O .5j} 

where Po is a constant, N s is the number of the scanning 
polygon sides, and h p is the pitch of the scanning locus 
7 e he,lx to to hel >cal reconstruction becomes a sequence 
of connected short lines which segments a helix (Fig 1 1) 
*k i 2 ShOWS four «»*nictBd images obtained using 

me three-dimensional broken-line cone-beam reconstruction 
algorithm, with po = 3. Ns = 8 and = 1.25 respectively 




E*JS- ^^nstrocted slices obtained using the three-dimensional 
broten-hne cone-beam reconstroction algorithm, when *> = 3. N s = 8. 
shown = ■ 0 ^j reSpeCtive,y ' COTTCS P° ndin 8 to the four slices of the phantom 



The averaged absolute errors between the original and recon- 
structed slices are 4.1 (upper left), 3.5 (upper right), 6.7 (lower 
left), and 5.7 (lower right) gray levels respectively. Similar 
remarks can be made as on the helical algorithm. The accuracy 
of the broken-line algorithm is very close to that of the helical 
algorithm. TTus algorithm will eventually be implemented for 
our x-ray microscope system. 

G. 3-D Dashed-Une Cone-Beam Reconstruction Algorithm 

The three-dimensional dashed-line algorithm is similar to 
the broken-line algorithm, except that the scanning turn of 
the three-dimensional dashed-line cone-beam algorithm is de- 
scribed as 



{ 



p = fp , 



(20) 



where p Q is a constant, N s is the number of the scanning 
polygon sides, and Ap is the counterpart of the pitches used in 
helical and broken-line algorithms. The difference between the 
broken-line algorithm and the dashed-line algorithm lies in that 
the inclined lines in the former become horizontal lines in the 
latter. Hence the scanning locus required by the dashed-line 
algorithm is similar to ''winding stairs" as shown in Fig. 13. 

The dashed-line algorithm is superior to the broken-line 
algorithm when the relative motion of a specimen along the z- 
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Fig. 13. The 



and a mid-plane. 




Fig. 14. Four reconstructed slices obtained using the three-dimensional 
dashed-line cone-beam reconstruction algorithm, when po = 3. Ns = &, 
and h p — 1.25. respectively, corresponding to the four slices of the phantom 
shown previously. 



axis is implemented by shifting the mechanical stage, allowing 
a much larger dynamic range for translating a rod-shaped spec- 
imen along the z -axis but producing greater mechanical errors. 
The dashed-line algorithm can greatly reduce the number of 
the mechanical stepping of the stage, and thus simplify the 
recalibration process needed for the reconstruction. 

Figure 14 shows four reconstructed images obtained using 
the three-dimensional dashed-line cone-beam reconstruction 
algorithm, with po = 3, Ns = 8 and h p — 1.25 respectively. 
The averaged absolute errors between the original and the 
reconstructed ones are 4.1 (upper left), 3.3 (upper right), 6.4 
(lower left) and 5.8 (lower right) gray levels respectively. 
Similar comments can be made as on the broken-line algo- 
rithm. This algorithm will also be implemented for our x-ray 
microscope system. 

H. Random Cone-Beam Reconstruction Algorithm 

We can also reconstruct the three-dimensional image of a 
specimen with a random scanning locus, which is modeled by 



P = Po + Cp * randomQ 
h = Ch • randomQ 



(21) 



algorithm with 
to the four 



where po. c p and are constants, and randomQ is a 
[-0.5,0.5] uniformly distributed random variable. Unlike in 
the previous cases where P was sampled evenly, 0 was 



Fig. 15. Four reconstructed slices obtained using the 
po = 3, c p as 1.0, and c h = 0.5, respectively, 
slices of the phantom shown previously. 



sampled randomly here. The random algorithm was designed 
for testing the stability of our cone-beam algorithm in case 
that p(P) does not describe a curve symmetric with respect 
to the origin. 

Figure 15 shows four reconstructed images obtained using 
the random algorithm with po = 3, c p = 1.0, and c h = 
0.5, respectively. The averaged absolute errors between the 
original and reconstructed slices are 4.0 (upper left), 13.4 
(upper right), 13.5 (lower left), and 13.4 (lower right) gray 
levels, respectively. Compared with the results obtained for 
Feldkamp's algorithm, the errors increased only slightly. This 
indicates that our cone-beam algorithm is quite robust. 



V. Analytic Properties of 
General Cone-Beam Algorithm 

Since our cone-beam algorithm is derived heuristically, like 
Feldkamp's algorithm, it is desirable to have some analytic 
properties. Assume that the reconstruction is performed using 
projection data collected from a turn of a scanning locus, three 
analytic properties of our algorithm are given in this section. 
These properties further justify the algorithm. 

A. Reconstruction Is Exact on Midplane 

When our cone-beam algorithm is applied in a planar 
scanning mode, a derivative-free noncircular fan-beam recon- 
struction formula is obtained by only considering the midplane 
reconstruction. It has been proved that the extended fan-beam 
formula is exact if some conditions are satisfied [41 J 1 . These 
conditions require that the scanning locus be differentiable 
almost everywhere with respect to the rotation angle, be 
symmetric with respect to the origin of the reconstruction 

1 An anonymous reviewer helped us correct an error contained in (23) in 
[41]. The work reported in [41] will be published in the October 1993 issue 
of IEEE Trans* on Image Processing. 
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coordinate system, and the mapping from fan-beam data to 
parallel-beam data is one-to-one (that is, the direction tangent 
to the locos stay outside of the fan delimited by the object 
support). At least for cone-beam x-ray micro tomography, these 
conditions are sufficiently general For example, polygons of 
even sides meet these requirements. 

The main idea used in the proof of [41] is that the extended 
fan-beam reconstruction formula can be expressed in terms 
of parallel-beam projection data. If the correct parallel-beam 
reconstruction formula could be obtained by rewriting the new 
fan-beam formula through transforming fan-beam data into 
parallel-beam data, then the new formula is proved 

The new formula is simpler than other published ones, as 
it is the same as the conventional equispatial one [12), except 
that the source-to-origin distance depends on the rotation 
angle. The new formula is also more stable, as it involves 
no derivative of the locus. 



nience, our cone-beam reconstruction formula is repeated here: 
2 Jo (P(0) - s) 2 

(24) 

We can rewrite the above formula as 



9(x,y, z) 



2 Jo (p(0)-s) 3 



f dpd0 (25) 



B. Vertical Integral of Reconstruction Is Exact where 

It has been proved that the vertical integral of the PSF 
sOr,y,*;ar 0 ,lto,zo) of our cone-beam reconstruction algo- 
rithm is a two-dimensional 6 function 6{x - x 0 )$(y - y Q ) in 
[42), where some intermediate integrals should be interpreted 
in the generalized sense when necessary [43], [12], [35], [10]. 
Using this interesting property of the PSF, it can be shown as 
follows that the vertical integral of a reconstructed 3-D image Cm the other hand, 
is equal to that of the real specimen image [42]. 

The reconstructed three-dimensional image can be expressed 
in terms of the spatially varying PSF as 



C=^» and *(/*) = *- 



aoJ J *( u,v » w M*»v»*; u > v > w ) d ^ vd «' 

(22) 

where i{x,y,z) is the real specimen image. Therefore, 
/ 9(x i y,z)dz 

J-oo 

/OO mOO fOO 
/ / »(u,v,w) 
oo J — oo J — oo 

/oo ^oo -oo 
oo /-oo J- ^ ^ ~ ~ V ) diidvdw 

~ I i{x y y,w)dw 

J —CO 

= / i{x,y y z)dz. 

J—oo 



t(ar, y, z) 

*Jq (p(/?)-s) 2 

J-oo 



-p) 



dpdp 



(23) 



C. Reconstruction Is Exact Ifi(x y y, z) Is Independent of z 

Let us assume that the two-dimensional scanning curve 
described by p(0),/3 € [0,2*] meet the requirements of our 
extended fan-beam reconstruction formula [41 J. For conve- 



(26) 

where T 0 (p t z) is the fan-beam projection data of the specimen 
shoe at the vertical axis coordinate z. S 0 (p, z) can be consid- 
ered as an approximation to Tpfaz). Spfaz) i s found by 
multiplying RpfaO with a projection factor V^Wh* ^ 
which is the cosine of the angle (depending on p) between fan 
rays corresponding to R 0 (p t f) and 7>(p, z) respectively. Note 
that i(x } y,z) is the true specimen image function, because 
our extended fan-beam reconstruction algorithm is exact in 
this case. When i(x,y>z) is independent of z t it is evident 
that Spfaz) = T 0 (p J z) t hence, the reconstruction is exact 
[44]. The above analysis also provides clues on why the helix- 
like scanning loci work better than multiple planar scanning 
loci. Generally speaking, Sp(p, z) obtained with a helix-like 
scanning locus approximates 7>(p,z) better. 

After a source point and a reconstruction point are given, 
there are infinite many ways to select tilted fans. The rationale 
behind our selection (and Feldkamp's selection as well) can be 
explained as follows. In our scheme, fans are tilted horizontally 
so that modified projection data can be imagined as being 
approximate to the imaginary horizontal fan-beam projection 
data, and if i(x,y t z) is independent of z y the reconstruction 
will be exact 
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VI. Conclusion and Discussion 

Our generalized cone-beam reconstruction algorithm pre- 
sented in this paper takes full advantage of the hardware 
characteristics of the x-ray microscope system being developed 
at SUNY at Buffalo, overcomes the limitations of Feld- 
kamp's cone-beam reconstruction algorithm and maintains its 
high computational efficiency and parallelism. Our cone-beam 
algorithm provides a framework allowing various scanning 
loci suitable for different applications. This algorithm not 
only gives satisfactory reconstruction images but also has the 
following three nice properties if a scanning locus meets some 
conditions. First, the reconstruction with a planar scanning 
locus gives the exact result on midplane. Second, 

g{x 1 y,z)dz- I %(x,y 1 z)dz (27) 

-oo J —oo 

where g(x t y,z) and i(x,y,z) are three-dimensional recon- 
structed and actual images respectively. Third, if i*(x, y, z) is 
independent of z, then 

g{x y y,z)^i(x,y,z). (28) 

Also, if a helix-like scanning locus is used, z-axis artifacts 
will be uniformized and reduced. 

Scanning loci used by the polygonal, multiple, dashed-line, 
and random algorithms do not satisfy Tuy's or Smith 1 s condi- 
tions for exact reconstruction [14], [351, [36], [10], therefore, 
an approximation cannot be avoided. Our approximation has 
been justified in this paper. On the other hand, helical and 
broken-line scanning loci do meet the requirements for the 
exact reconstruction. Our algorithm does not take advantage 
of this theoretical potential. The exact reconstruction schemes 
[35], [36], [10] could be used in this case. However, these 
exact algorithms are computationally inefficient. On the other 
hand, as shown in this paper, our approximate algorithms are 
quite accurate. Although Smith also proposed a both exact and 
efficient scheme [10, 33], as mentioned earlier, the conditions 
for his exact and efficient algorithm are too restrictive. 

Our cone-beam reconstruction algorithm can be considered 
as being derived based on our derivative-free noncircular fan- 
beam reconstruction formula [41]. Using our methodology, 
two other noncircular fan-beam reconstruction formulas [28], 
[29], [45], [46] can also be extended to cone-beam recon- 
struction with various non-planar scanning loci. Although a 
derivative of the scanning locus will be contained in these 
cone-beam reconstruction formulas, it can be readily and 
similarly proved that all three analytic properties described 
in Section V still hold. By doing so, we actually extend Gull- 
berg's cone-beam reconstruction formula [29] into a general 
case. In fact, Gullberg et aL extended Feldkamp's algorithm 
horizontally* while we did it both horizontally and vertically. 

In order to study the reconstruction accuracy of our cone- 
beam algorithm, a systematic and quantitative error analysis is 
being performed [44], [47]. Relevant results will be reported 
in the future. 
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